clear;
n = 100;
t = 0;
dt = 0.001;
nt = 10;
beta = 4;
dx = 1./n;
qa = zeros(n,1);
q = qa;
f = zeros(n+1,1);
x = q;

for i = 1:n
  x(i) = (i-0.5)*dx;
  q(i) = 4;
  if (i>n/2) q(i) = 2; endif
endfor

plot(x,q,"-"); pause

for m = 1:nt
  # autalizamos
  t = t+dt;
  qa = q;
  
  # calculo dos fluxos
  for i=2:n
    # Option 1: centrado
    #qq = (qa(i-1)+qa(i))/2;
    
    # Option 2: full upwind, Godunov
    if (beta>0) 
      qq = qa(i-1);
    else qq = qa(i); endif;
    
    f(i) = beta*qq;
  endfor

  f(1) = beta*qa(1);
  f(n+1) = beta*qa(n);

  # calculo do q
  for i=1:n
    q(i) = qa(i)-dt/dx*(f(i+1)-f(i));
  endfor
  
    
  # graficamos
  plot(x,q,"@"); pause;
endfor

